# With Village Fixed Effects ---------------------------------------------------

# LPG Connection (Binary)
model.dd.vfe.uselpg.bplaay = fixest::feols(uselpg ~ bplaay*wave + age + gender +
                                           religion + education + hhsize +
                                           lnexp| village, 
                               data = pool, weights = pool$weights)

# LPG Home Delivery (Binary)
model.dd.vfe.convenience.bplaay = fixest::feols(convenience ~ bplaay*wave + age + gender +
                                           religion + education + hhsize +
                                           lnexp  | village, 
                                         data = pool, weights = pool$weights)

# Grid Connection (Binary)
model.dd.vfe.usegrid.bplaay = fixest::feols(usegrid ~ bplaay*wave + age + gender +
                                            religion + education + hhsize +
                                            lnexp  | village, 
                                          data = pool, weights = pool$weights)

# With Household Fixed Effects -------------------------------------------------

# Regression Models for Energy Access Response Variables
# Robust standard errors are clustered at household level.

# LPG Connection (Binary)
model.dd.hfe.uselpg.bplaay = fixest::feols(uselpg ~ bplaay*wave + hhsize +
                                           lnexp  | hh, 
                                         data = pool, weights = pool$weights)

# LPG Home Delivery (Binary)
model.dd.hfe.convenience.bplaay = fixest::feols(convenience ~ bplaay*wave + hhsize +
                                                lnexp  | hh, 
                                              data = pool, weights = pool$weights)

# Grid Connection (Binary)
model.dd.hfe.usegrid.bplaay = fixest::feols(usegrid ~ bplaay*wave +  hhsize +
                                            lnexp  | hh, 
                                          data = pool, weights = pool$weights)

## Model and coefficient lists for village and household fixed effects combined ##
models.vfe.hfe.dd_access_bplaay <- list(model.dd.vfe.uselpg.bplaay, 
                                  model.dd.vfe.convenience.bplaay,
                                  model.dd.vfe.usegrid.bplaay, model.dd.hfe.uselpg.bplaay,
                                  model.dd.hfe.convenience.bplaay, 
                                  model.dd.hfe.usegrid.bplaay)

esttex(models.vfe.hfe.dd_access_bplaay, 
       se = "cluster", 
       digits = 3,
       dict = c(bplaay = "BPL/AAY", wave = "2018", uselpg = "LPG Connection",
                usegrid = "Grid Connection", convenience = "LPG Home Delivery",
                village = "Village", hh = "Household"),
       keep = c(":", "2018"),
       fitstat = c("r2"),
       fixef_sizes = T,
       float = F,
       replace = T,
       file = c("./Manuscript/Tables/dd_access_bplaay.tex"))